from libc.stdio cimport *

cimport cython
import numpy as np
cimport numpy as np
 
cdef extern from "stdio.h":
    FILE *fopen(const char *path, const char *mode)
    int fclose(FILE *stream)
    ssize_t getline(char **lineptr, size_t *n, FILE *stream)

cdef extern from "string.h":
    char *strchrnul(const char *s, int c)

cdef extern from "stdlib.h":
    double atof(const char *nptr)
    void *malloc(size_t size)
    void free(void *ptr)

@cython.boundscheck(False)  # Deactivate bounds checking
@cython.wraparound(False)   # Deactivate negative indexing.
def covMatrix(filename,int ndim):
    """
	C优化IO和numpy.outer
	"""
    filename_byte_string = filename.encode("UTF-8")
    cdef char* fname = filename_byte_string
 
    cdef FILE* cfile
    cfile = fopen(fname, "rt")
    if cfile == NULL:
        raise FileNotFoundError(2, "No such file or directory: '%s'" % filename)
 
    cdef char * buf = <char *>malloc(1024*sizeof(char))
    cdef size_t maxlineSize = 1024
    cdef char * p = buf
    cdef ssize_t read
 
    r = np.zeros(ndim,dtype=np.double)
    cdef double [:] r_view = r
    
    cdef Py_ssize_t  i = 0
    cdef int n = 0
    s1 = np.zeros(ndim,dtype=np.double)
    cdef double [:] s1_view = s1
    s2 = np.zeros((ndim,ndim),dtype=np.double)
    cdef double [:,::1] s2_view = s2
  
    while True:
        read = getline(&buf, &maxlineSize, cfile)
        if (read <= 0): break
        p = buf
        for i in range(ndim):
            r_view[i] = atof(p)
            p = strchrnul(p,',')        
            p +=1
        n += 1 
        for i in range(ndim):
            s1_view[i] += r_view[i]
            for j in range(i,ndim):
                s2_view[i][j] += r_view[i]*r_view[j]   #保留上三角
#        s1 += r
#        s2 += np.outer(r, r)
    free(buf)    
    fclose(cfile)

    s1 /= n
    s2 /= n

    for i in range(ndim):
        for j in range(i,ndim):
            s2_view[i][j] -= s1_view[i]*s1_view[j]  #保留上三角
#    s2 -= np.outer(s1, s1) #总体协方差

    #s2 *= (n/(n-1))   #样本协方差，无偏估计
    d = np.sqrt(s2.diagonal())
    cdef double [:] d_view = d    

    for i in range(ndim):
        for j in range(i,ndim):
            s2_view[i][j] /= d_view[i]*d_view[j]
            s2_view[j][i] = s2_view[i][j]           #补充下三角
#    s2 /= np.outer(d, d)  #归一化相关系数矩阵
 
    return s2